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Abstract 

In  this  report,  we  discuss  and  compare  several  methods  for  polynomial  Interpolation  of 
Global  Positioning  System  ephemeris  data. 

Intoduction 

The  problem  of  interpolating  Global  Positioning  System  (GPS)  ephemeris  data  is  an 
important  aspect  of  GPS.  Given  that  a  high  accuracy  (<  Im  ),  high  precision  (1  cm)  orbit 
can  be  generated  though  the  use  of  dense  observations  and  special  integrations,  it  is  necessary 
to  interpolate  these  ephemeris  at  high  accuracy  to  utilize  these  orbits. 

These  high  accuracy  orbits  are  produced  by  several  organizations  (DMA,  NGS,  JPL, 
several  Universities)  and  are  widely  available.  An  ephemeris  typically  consists  of  satellite 
positions  at  evenly  spaced  times  over  a  week.  Most  ephemeris  are  given  at  900  sec  (15  min) 
time  steps  although  the  NGS  ones  are  at  1200  sec  (20  min).  The  GPS  satellites  are  in  12  hr 
circular  orbits  making  900  sec  ephemeris  steps  7.5  deg  of  arc. 

The  typical  user  collects  GPS  data  at  intervals  from  1  sec  to  30  sec  and  needs  to  find 
the  satellite  position  at  the  times  of  that  data.  The  times  needed  are  really  not  the  evenly 
spaced  received  times,  but  the  transmit  times  that  are  about  60  msec  before  reception.  A 
precise  value  for  this  propagation  delay  is  not  known  until  the  solution  process  in  partially 
done.  Therefore  usually  one  needs  to  find  a  cluster  of  satellite  positions  a  few  msec  from  a 
nominal  evenly  spaced  interval. 

In  the  past  the  typical  technique  used  by  the  DMA  [Malys,  1989],  NGS  [Remondi,  1991], 
JPL  [Watkins,  1995]  and  others  is  a  Lagrange  interpolation.  The  orders  vary  from  to 
11*^.  This  approach  directly  computes  the  value  of  the  function  (the  three  Cartesian  Earth 
centered  earth  fixed  coordinates)  from  the  unique  polynomial  going  through  the  data  points. 
The  coefficients  are  not  found,  and  finding  them  may  introduce  errors  [Press  et  al,  1992]. 
Several  evaluations  of  the  accuracy  of  this  method  [Remondi,  1991,  Smith  and  Curtis,  1983] 
have  been  made.  It  is  generally  found  that  an  8*^  order  Lagrangian  interpolation  using 
900  sec  data  with  the  unknown  in  the  center  of  the  points  gives  values  that  compare  with 
numerical  integration  at  the  1  cm  level. 

The  problem  addressed  here  is  to  find  if  a  more  efficient  numerical  method  that  achieves 
the  same  accuracy  can  be  used.  This  is  motivated  by  the  movement  of  processing  from 
mainframes  to  PC’s  (486’s  and  above). 

Several  aspects  unique  to  the  GPS  satellites  make  this  problem  of  interpolating  the 
data  different  from  the  general  problem  of  interpolation.  Though  we  are  interpolating  GPS 
satellite  orbital  position  data,  it  may  be  that  the  methods  here  are  applicable  to  a  broader 
class  of  problems.  Where  possible,  we  intend  to  take  full  advantage  of  the  special  geometry 
of  the  GPS  satellite  orbits. 
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Figure  (1)  -  4  day  x-coordinate  GPS  ephemeris  data  at  900  second  intervals 

A  typical  precise  ephemeris  orbit  is  supplied  over  an  interval  of  eight  days.  Each  ephemeris 
overlaps  1  day  at  each  end  with  another  ephemeris.  It  consists  of  position  data  which  is 
Earth-centered,  Earth-fixed  Cartesian  coordinates  given  every  900  seconds  (t,x,y,z).  We 
will  not  include  velocity  data  which  may  in  some  cases  be  available.  A  plot  of  the  data  shows 
that  it  is  “almost”  periodic  with  a  period  of  24  hours  or  96  intervals  of  900-second  each;  Of 
course  the  data  would  be  periodic  in  inertial  coordinates.  The  data  is  given  in  a  rotating 
frame  (Earth  centered  -  Earth  fixed)  however.  This  is  done  to  place  several  subtle  effects, 
such  as  polar  motion,  into  the  ephemeris  generation  problem.  The  user  then  does  not  have 
to  have  access  to  current  data  for  geophysical  effects  to  compute  an  Earth  fixed  position 
.solution.  Figure  (1)  is  a  plot  of  x-coordinate  data  (in  kilometers)  over  a  four  day  period. 
Note  that  the  plot  is  with  respect  to  the  node  (point)  number  which  ranges  from  1  to  384.  It 
is  convenient  to  use  node  numbers  since  we  may  choose  to  map  the  interval  of  interest  into 
different  subintervals.  For  instance,  the  Lagrange  interpolation  method  maps  the  interval 
of  interest  into  the  interval  [—1,1)  while  the  trigonometric  polynomial  interpolation  method 
maps  the  interval  of  interest  into  the  interval  [0, 27r).  Later  we  will  use  the  node  numbers  to 
compare  the  residuals  for  different  methods. 

There  are  two  lengths  of  computation  time  in  which  we  are  interested  for  this  problem. 
One  time  length  is  the  generation  time  of  the  function  (interpolant)  which  interpolates  the 
evenly-spaced  data.  The  other  is  the  time  needed  to  evaluate  the  interpolant  at  a  point.  For 
polynomial  interpolation  these  might  be  the  calculation  of  the  coefficients  of  the  polynomial 
interpolating  the  data  and  the  time  needed  to  evaluate  that  polynomial,  once  generated,  at 
a  particular  point.  As  we  shall  see,  however,  there  are  clever  ways  in  which  to  minimize  the 
amount  of  work  done  in  generating  the  desired  data,  and  the  two  times  might  not  be  easily 
distinguished  in  these  cases.  In  recommending  a  particular  technique  for  interpolation  it  is 
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important  to  know  whether  the  interpolant  will  be  evaluated  once  (or  just  a  few  times)  on  an 
interval  or  if  it  will  be  evaluated  many  times  throughout  the  interval.  With  a  cluster  of  times 
to  be  interpolated  on  an  interval,  the  cost  of  generating  the  interpolant  should  be  amortized 
over  the  set  of  times.  This  means  that  the  time  needed  to  generate  the  interpolant  (a  one¬ 
time  cost)  would  not  be  as  much  a  concern  as  the  time  needed  to  evaluate  the  interpolant 
at  each  of  the  desired  times  (a  recurring  expense).  On  the  other  hand,  if  one  or  a  very  few 
points  were  to  be  calculated  on  a  given  interval,  the  time  needed  to  generate  the  interpolant 
would  probably  be  more  significant  a  concern  than  the  time  needed  to  evaluate  it  at  the 
desired  times  since  in  general  the  generation  of  an  interpolant  is  much  costlier  (in  time)  than 
its  evaluation.  This  is  a  similar  argument  as  one  finds  in  deciding  whether  to  use  Gaussian 
elimination  or  LU  factorization  in  solving  systems  of  linear  equations. 

One  method  of  polynomial  interpolation  involves  performing  the  interpolation  at  a  point 
without  actually  calculating  the  coefficients  of  the  interpolating  polynomial.  This  involves 
many  less  operations  than  the  evaluation  of  a  polynomial  of  degree  n  at  a  particular  point. 
This  is  a  common  technique  currently  employed.  If  we  do  not  explicitly  calculate  the  in¬ 
terpolant,  however,  we  will  in  general  still  need  to  calculate  some  quantities  prior  to  in¬ 
terpolation.  For  instance,  we  shall  see  that  divided  differences  would  need  to  be  available 
prior  to  using  the  nested  multiplication  algorithm  used  in  evaluating  the  Newton  form  of  the 
interpolating  polynomial. 

Since  ephemeris  data  is  generated  for  an  eight  day  time  period,  we  have  the  opportunity 
to  “front  load”  our  work  at  the  time  of  ephemeris  receipt.  By  calculating  needed  data 
in  advance  we  should  be  able  to  shorten  the  real  time  operation  count.  Thus  we  will  be 
more  concerned  with  the  rapid  evaluation  of  interpolants  for  specific  times  than  their  rapid 
generation.  Of  course  in  some  cases  the  times  of  evaluation  and  generation  will  be  closely 
related,  as  mentioned  above.  In  others,  they  will  be  quite  different  and  our  hope  will  be  to 
shift  as  much  of  the  work  as  possible  to  the  generation  of  efficient  interpolants  so  as  to  allow 
the  rapid  evaluation  of  those  interpolants. 

In  summary,  we  will  describe  the  following  methods: 

•  Lagrange  polynomial  interpolation 

•  Newton’s  divided  difference  interpolation  polynomial 

•  Difference  Tables 

•  Cubic  Spline  interpolation 

•  Trigonometric  polynomial  interpolation 

•  Tshebyshev  polynomial  interpolation 

We  will  describe  the  advantages  and  disadvantages  of  each  of  these  methods  for  the  problem 
of  interest,  namely  for  the  efficient  interpolation  at  clusters  of  points. 
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Lagrange  Polynomial  Interpolation 

Before  we  begin  our  investigation,  it  is  necessary  to  describe  the  method  which  is  currently 
being  used.  Simply  put,  given  the  n  + 1  ephemeris  values  /(to),  •  •  •  ,f{tn)  at  the  distinct 

times  to,  ti, . . . ,  t„,  there  exists  a  unique  interpolating  polynomial  pn  satisfying 

Pn(tt)  =  / (tj),  Z  =  0,  1,  .  .  .  ,  U 

This  polynomial  can  be  written  in  the  form  (called  the  Lagrange  interpolation  polynomial) 

Pn{t)  =  t.mw  (') 

t=0 


m  = 


where 

(t  —  to)(t  —  ti)  •  •  •  (t  —  t,-i)(t  —  ti+i)  •  •  •  (t  —  t„) 

(tj  —  to)(tj  ti)  •  •  •  (tj  t,_i)(tj  tj-f-l)  •  •  •  (tj  tn) 

The  eleventh  order  Lagrange  method  uses  twelve  data  points  to  generate  an  eleventh  order 
polynomial  according  to  equation  (1).  This  polynomial  can  then  be  evaluated  at  desired  times 
within  the  interval  of  interest.  The  error  Rn{t)  in  using  the  Lagrange  interpolant  p„(t)  to 
estimate  the  function  /(t)  (having  at  least  n  +  1  derivatives  throughout  the  open  interval) 
at  some  point  t  can  be  written  [Buchanan  and  Turner,  1992]: 


Rn{t)  -  f{t)  -  Pnit)  = 

where  if  is  some  point  in  the  interval  [to,  tn]  and 


-'n+l 


(t)/("+i)(f) 


(u  +  l)! 


i-0 

One  difficulty  in  implementing  high  degree  polynomial  interpolation  routines  of  any  kind 
is  the  fact  that  the  error  between  the  interpolating  polynomial  and  the  data  or  function  being 
interpolated  grows  rapidly  near  the  endpoints  of  the  interval  over  which  the  interpolation 
is  being  performed.  For  this  reason  the  eleventh  order  Lagrange  method  is  overlapped  as 
successive  intervals  are  chosen  within  the  ephemeris  (we  call  this  walk-along  interpolation). 
Due  to  the  high  accuracy  requirements,  only  the  center  subinterval  is  interpolated  for  each 
Lagrange  polynomial  which  is  generated.  Whereas  the  initial  interval  spans  points  one 
through  twelve,  the  second  interval  spans  points  two  through  thirteen  in  order  to  provide 
the  highest  degree  of  accuracy.  The  first  eleventh  order  Lagrange  polynomial  would  then 
be  used  for  times  between  points  six  and  seven,  while  the  second  polynomial  would  be  valid 
for  times  between  points  seven  and  eight.  The  numerical  accuracy  of  this  method  has  been 
verified  to  the  1  cm  level  for  the  data  we  are  interpolating  [Remondi,  1991). 

Difficulties  arise  in  that  the  process  of  creating  and  evaluating  the  resulting  eleventh 
order  polynomials  is  computationally  expensive.  The  cost  of  evaluating  the  Lagrange  form 
at  a  point  is  provided  by  de  Boor  [1978].  It  is 


{2n-2)A  +  {n-2)M  +  {n-l)D 
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for  each  of  the  n  +  1  numbers  where  A  denotes  an  addition  or  subtraction  and  M  and 
D  denote  multiplication  and  division,  respectively.  Forming  equation  (1)  then  takes  another 

(n  —  1)  A  +  n  M 

operations,  leaving  the  total  count  per  component  of  the  position  vector  at 

(3n  -3)A  +  (2n-2)M +  (n- 1)D 


This  is  the  number  of  operations  per  component  in  the  implementation  of  Lagrange  inter¬ 
polation.  A  simple  modification  of  the  algorithm  would  reduce  the  amount  of  work  to 

(2n  —  l)A-hnMAnD 

operations  (see  de  Boor  [1978]).  It  consists  of  first  forming  the  quantity 

m 


Ilj^i  ^j) 

Afterwards,  Pn{t)  is  calculated  through 


z  =  1 , . . . ,  n 


n 

=  n 

t=l 

Pn{i)  =  4>n{t)  ^  Yi/{t  -  ti) 
i=l 

This  method  is  somewhat  faster  than  the  method  currently  being  used  and  would  be  easy 
to  implement  (see  program  lagrange  in  Appendix  A).  In  addition,  no  loss  of  accuracy  would 
occur  in  its  implementation.  If  the  time  t  is  very  close  to  one  of  the  interpolating  points  t,, 
one  must  be  careful  in  computing  Pn{t)  because  of  the  division  of  by  a  very  small  number. 


Newton’s  Divided  Difference  Interpolation  Polynomial 


A  more  efficient  means  by  which  we  may  form  the  interpolating  polynomial  is  through 
divided  differences.  We  follow  the  developments  given  in  Buchanan  and  Turner  [1992]  and 
de  Boor  [1978].  Suppose  we  have  n  -f- 1  distinct  interpolation  points  to,  ti, . . . ,  Define  the 
zeroth  divided  difference  at  t,  by  /[f,]  =  /(L).  The  first  divided  difference  at  ti,tj  is  defined 
by 


(Xq  —  fli'i, 


/[L]  -  fitj] 

ti  tj 


and  the  kth  divided  difference  f[to,ti, . . . ,  tfc]  is  defined  by 


flfc  —  /[^Oj  fl?  •  •  •  5  tk\  — 


/[to,  fl,  .  .  ■  ,  tfc-l]  -  f[tl,t2,  ...,tk\ 
i-o  —  i-k 


fc>0 
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Newton’s  divided  difference  interpolation  polynomial  is  the  interpolation  polynomial  agreeing 
with  the  function  /  at  the  points  to,t\, . . .  ,tn  is  given  by 

Pn{t)  =  ao  +  (t  —  ^o)fll  +  (^  ~  io){l  ~  ^l)®2  +  •  •  •  +  (^  ~  ^o)(^  —  t\)  •  ■  ■  {t  —  tn-l)an  (2) 
or,  rearranging, 

71-1  times 

=  cto  "h  ~  ^o){®l  +  (^  ~  ^l){^2  +  •  ■  ■  +  ~  tn-2}{^n-l  +  ~  ^n-l)^n  }•••}} 

This  form  consists  of  two  additions  and  a  multiplication  per  level  in  the  expression.  Since 
there  are  n  levels,  the  operation  count  is  seen  to  be 

n{2A  +  M) 

which  is  more  economical  than  the  standard  Lagrange  form.  This  leads  us  to  the  so-called 
nested  multiplication  or  Horner’s  algorithm: 

Given  the  n  -|-  1  distinct  points  to^ti, . , .  ffn  with  associated  coefficients  ao,  ui, .  •  ■ ,  Un,  the 
value  of  the  interpolating  polynomial  p„(t)  for  some  t  G  is  given  by  bo  according  to 

the  following  iteration: 

Set  bji  —  O/ji 

For  k  =  n  —  I  to  0  by  -1 

bk  =  o-k  +  {t  ~  bk+i 

End  For 

By  the  uniqueness  of  the  interpolating  polynomial  there  can  be  no  difference  in  compar¬ 
ison  to  Lagrange,  but  the  gain  in  speed  may  be  of  importance  for  our  purposes.  Note  that 
the  divided  differences  Ufc  can  be  calculated  and  stored  in  advance  of  the  actual  interpola¬ 
tion  so  that  the  operation  counts  given  here  reflect  the  operations  needed  at  interpolation 
time.  A  total  operation  count  would  have  to  include  the  operations  needed  to  generate  the 
divided  differences.  Also,  calls  from  storage  may  need  to  be  counted,  depending  on  system 
architecture  considerations. 

The  divided  difference  algorithm  does  not  take  advantage  of  the  fact  that  our  interval 
sizes  are  fixed  —  we  required  the  nodes  to  be  distinct  but  made  no  restriction  on  the  spacing 
between  nodes.  In  the  next  section  we  will  investigate  the  special  case  when  the  interval 
sizes  are  constant. 


Difference  Tables 

The  case  of  equally-spaced  data  points  is  a  special  case  of  Newton’s  divided  differences 
and  leads  to  other  interpolation  formulas.  The  error  and  operation  counts  for  the  methods 
presented  here  are  essentially  identical  to  those  presented  above.  The  formulas  are  given 
in  their  simplest  form  and  should  not  be  used  for  computation.  A  nested  multiplication 
approach  similar  to  the  one  described  in  the  previous  section  should  be  used  for  each  of 
these  in  order  to  minimize  the  cost  of  computation.  One  important  aspect  of  this  method 
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is  the  determination  of  the  differences  and  the  method  to  be  used  to  interpolate  at  a  given 
time  t.  It  will  be  necessary  to  include  some  code  to  determine  which  differences  are  to  be 
used,  though  the  differences  themselves  can  be  calculated  when  the  given  ephemeris  becomes 
available.  In  addition,  the  chosen  method  will  depend  on  the  location  of  the  interpolation 
time  relative  to  the  data  times.  Here  we  follow  the  description  given  in  Buchanan  and  Turner 
[1992].  Our  data  occurs  at  times  which  can  be  expressed  as 

4  =  *0  +  kh 

where  to  is  a  reference  time  for  the  interval  of  interest  and  h  is  the  constant  steplength.  We 
normally  think  oik  zs  being  a  positive  integer  and  to  as  being  the  initial  time  in  the  interval 
of  interest  but  in  this  case  we  will  only  require  k  to  be  an  integer  and  to  to  be  any  time 
corresponding  to  a  data  point  in  the  interval.  The  sign  of  k  will  depend  on  the  reference 
time  to  in  relation  to  the  time  of  interest.  There  are  now  several  differences  which  can  be 
defined,  one  of  which  is  the  forward  difference. 

The  general  forward  difference  Af{ti)  is  given  by 

A/(ti)  =  f{ti+i)  -  f{ti)  =  fiU  +  h)-  f{ti) 

Its  powers  are  calculated  recursively  according  to  the  following 

A‘/((,)  =  A(A‘-'/(t,))  =  A‘-'/(ii+i)  -  A‘-V(e,-) 


j=0  \J/ 


where  we  have  introduced  the  notation  fi  =  Additionally,  the  differences  are  related 

to  divided  differences  by 

A>‘f{U)  =  k\h^f[tiffi+u---U+k] 

An  application  of  this  last  formula  to  equation  (2)  immediately  yields  a  forward  difference 
-formula  (called  Newton’s  forward  difference  formula  or  the  Newton-Gregory  forward  differ¬ 
ence  formula).  Here  we  assume  that  the  degree  of  the  interpolating  polynomial  is  n  while 
the  number  of  data  points  in  our  table  is  N: 


Pn{t)  =  /o  +  - ^-^A/o  + 


it-to){t-h) 


AVo  +  ---  +  ^ 


{t  -  to){t  -  ti)---{t-  f„_i) 


A  simple  change  of  variable  t  —  (t  —  to)/h  yields  the  compact  form 

Pnir)  = 

j=o  \d/ 

with  the  generalized  binomial  coefficients 


AVo 


r(r  -  1)  •  •  •  (r  -  j ■  +  1) 


Examination  of  these  formulas  reveals  that  their  operation  counts  could  be  made  similar 
to  the  operation  count  for  the  Newton  divided  difference  formula  with  a  constant  step  size 
through  the  use  of  nested  multiplication.  Their  accuracy  is  the  same  as  the  Newton  divided 
difference  formula  so  that  none  should  be  ruled  out  in  the  search  for  a  more  rapid  method 
of  interpolation. 


Cubic  Spline  Interpolation 


Cubic  spline  interpolation  is  computationally  efficient  and  has  an  advantage  with  respect 
to  walk-along  Lagrange  because  it  allows  the  user  to  calculate  the  interpolating  polynomials 
over  the  entire  interval  at  one  time,  at  the  beginning  of  the  interpolation  process.  This 
calculation  involves  solving  a  tridiagonal  system  of  equations.  Additionally,  the  fact  that  each 
subinterval  is  represented  by  a  cubic  polynomial  means  that  evaluations  on  those  intervals 
are  much  quicker  than  their  eleventh  order  polynomial  counterparts,  making  the  cubic  spline 
method  a  good  choice  where  accuracy  is  less  important  than  speed. 

Here  we  will  sketch  the  derivation  of  the  cubic  spline  interpolation;  thorough  treatments 
of  cubic  spline  interpolation  can  by  found  in  [Ahlberg,  Nilson  and  Walsh,  1967],  [Buchanan 
and  Turner,  1992],  [de  Boor,  1978],  and  [Press,  1992].  Given  the  n  +  1  ephemeris  values 

at  the  distinct  times  a  =  =  b,  we  construct  a  piecewise 

cubic  interpolant  p  to  /  as  follows.  On  each  subinterval  [<„  U+i]  we  wish  to  construct  a  cubic 
polynomial  p,-  in  such  a  way  that  the  resulting  interpolation  formula  over  the  entire  interval 
is  continuous  in  its  first  and  second  derivatives.  The  result  is 

p.(t)  =  Mt)  fiU)  +  Bi{t)  /(t.+i)  +  Ci{t)  f"{u)  +  Di{t)  r (i.+i),  i  =  0, . . . ,  n  -  1  (3) 


where 


A.  =  ^ 

’  ti+i-ti’ 


t  —  U 


We  do  not  yet  have  the  n  +  1  values  of  f"{U)  needed  for  the  determination  of  the  solution, 
but  application  of  the  continuity  of  the  first  derivative  on  the  entire  interval  leads  us  to  the 
following  equations: 


ti  —  ^  ii+l  ii-1  fU/j,  \  .  ^i+1  riif,  \  /(^«)  /(^») 

(ti-l)  +  - 5 - /  {ti)  +  {ti+l)  -  —7-1737;  t.  - 


U  — 


(4) 

for  z  =  1,2, ...  ,n  -  1.  Note  that  there  are  n  -  1  equations  for  n  +  1  unknowns,  leaving 
two  second  derivatives  undetermined.  The  choice  of  the  two  boundary  conditions  /"{a)  and 
f"{h)  provides  the  required  unique  solution.  In  our  case  it  makes  sense  to  use  the  periodicity, 
i.e.  apply  (3)  and  (4)  for  z  =  n  and  enforce  po(t)  =  Pn{t  +  ^n)  for  to  <  t  <  ti.  As  a 
consequence,  /(to)  =  f{tn),  f{ti)  =  /(Wi)?  f"{to)  =  f"{tn),  and  /"(ti)  =  /"(tn+i)- 
Therefore  if  we  use  equal  step  size  /i,  the  cubic  spline  version  of  equation  (4)  can  be  written 
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as 

/4  1  0  0 

14  10 

0  14  1 

0  •••  0  1 

\1  •••  0  0 

where  /?,■  =  24/[i,_i,i,,t,+i]. 

If  accuracy  can  be  sacrificed  for  speed  then  the  cubic  spline  method  may  be  preferred 
over  any  of  the  methods  here.  However,  the  0{h‘^)  accuracy  provided  by  the  cubic  spline 
may  be  insufficient  for  GPS  satellite  interpolation  requirements. 

Trigonometric  Polynomial  Interpolation 

The  roots  of  this  method  of  interpolation  can  be  traced  to  the  beginning  of  the  nineteenth 
century.  Briggs  and  Henson  [1995]  present  a  brief  history  of  this  method,  in  particular  the  fact 
that  Gauss  used  it  around  1800  to  interpolate  the  orbit  of  the  asteroid  Pallas.  The  preceding 
methods  are  standard  interpolation  techniques  typically  used  for  continuous,  differentiable 
functions  defined  on  compact  intervals.  No  other  special  information  about  the  functions 
being  interpolated  is  exploited  by  these  methods.  It  is  at  this  point  that  we  examine  some 
special  properties  of  our  GPS  ephemeris  data.  As  previously  mentioned,  our  ephemeris 
data  is  supplied  over  an  interval  of  eight  days  and  consists  of  Earth-centered,  Earth-fixed 
Cartesian  coordinate  position  data  given  every  900  seconds.  A  plot  of  the  data  in  Figure 
(1)  shows  that  it  is  very  close  to  periodic,  and  it  is  for  this  reason  that  we  examine  the 
trigonometric  polynomial  interpolation  method. 

Due  to  the  fact  that  the  position  data  has  a  period  of  twenty-four  hours,  we  restrict  our 
attention  to  a  single  twenty-four  hour  period  and  generate  a  trigonometric  polynomial 
using  all  the  data  available  over  that  period.  Again,  we  must  remember  that  our  satellite 
orbit  is  not  truly  periodic,  but  very  close  to  that  in  inertial  coordinates.  Since  the  error 
incurred  by  assuming  the  data  to  be  periodic  over  a  j  day  period  would  be  almost  j  times 
as  great  as  the  error  incurred  from  assuming  the  orbit  to  be  periodic  over  a  single  day,  we 
discourage  use  of  this  routine  over  intervals  exceeding  the  fundamental  period  of  the  data. 
In  order  to  minimize  the  effects  of  the  assumption  of  periodicity  one  should  interpolate  over 
a  single  period. 

The  idea  in  generating  our  interpolant  is  to  assume  that  the  data  is  from  a  periodic 
function  of  time  defined  over  the  interval  [0, 27r)  which  can  be  represented  by  a  trigonometric 
polynomial  of  the  form 


4 

1 


1 

4 


\ 

rih)  \ 

/  Ri  \ 

nt2) 

R2 

f"{h) 

Rz 

ntn-x) 

/ 

f'itn)  } 

Rn  J 

(5) 


n 

Pn{t)  =  ao  -f  ^  (ofc  cos  kt  +  hk  sin  kt) 
k=l 

In  our  problem  we  simply  map  the  twenty-four  hour  interval  of  interest  into  the  interval 
[0, 27r)  using  a  linear  change  of  variable.  In  deriving  the  coefficients  we  follow  the  discussion 
of  interpolating  polynomials  and  the  Fast  Fourier  Transform  found  in  Buchanan  and  Turner 
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[1992].  Clearly,  the  above  equation  can  be  written  as 

Pn(t)=  £ 

k=—n 

where  the  coefFcients  are  given  by 


7fc  =  «*  +  if^k  k  =  -ra, . . . ,  n 

If  we  denote  the  point  on  the  unit  circle  corresponding  to  t  €  [0, 27r)  by 

e  = 


we  may  then  write 

Pnit)=  ^ 
k——n 

If  we  denote  the  points  corresponding  to  nodes  by 

=  cos  tj  +  i  sin  tj 

the  interpolation  problem  is  to  find  the  coefficients  satisfying 

p.(ti)=  E  >{?  =  /(« 

k=—n 

where  f{tj)  are  the  function  values  for  j  =  0, 1, . . . ,  2n.  For  an  odd  number  of  nodes,  it  can 
be  shown  that 

1 

. " 

Unfortunately,  as  pointed  out  by  Buchanan  and  Turner  [1992],  the  operation  count  for 
this  discrete  Fourier  Transform  is  O(n^)  making  it  too  expensive  for  practical  application. 
Luckily,  there  is  a  faster  way  to  calculate  the  discrete  Fourier  coefficients. 

Suppose  that  we  have  the  case  where  the  number  of  nodes  is  of  the  form  n  =  2^ .  The 
task  is  then  to  find  the  trigonometric  polynomial  which  interpolates  /  at  the  uniformly 
spaced  data  set  In  tkis  case  the  trigonometric  polynomial  can  be  calculated 

by  using  the  Fast  Fourier  Transform,  which  has  an  operation  count  of  0(n  log  n).  This 
makes  it  an  attractive  method  for  our  purposes  arid  we  will  refer  to  it  as  the  FFT  method. 
Unfortunately,  however,  over  a  twenty-four  hour  period  our  data  will  consist  of  ninety-six 
points  separated  by  intervals  of  900  seconds.  Though  ninety-six  is  not  a  power  of  two,  we 
nevertheless  can  use  n  =  2®  =  32  points  in  the  interval  to  test  our  trigonometric  polynomials. 
We  simply  select  every  third  point  in  the  interval  and  calculate  the  trigonometric  polynomial 
which  interpolates  at  those  thirty-two  points.  We  may  then  compare  the  Lagrange  method 
to  the  FFT  method  by  choosing  some  subinterval  having  a  length  of  36  nodes,  twelve  of 
which  are  used  to  generate  the  eleventh  order  Lagrange  interpolating  polynomial.  These 
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twelve  nodes  (equally  spaced  with  3  times  the  spacing)  should  coincide  with  twelve  of  the 
thirty-two  nodes  used  in  the  FFT  method. 

The  above  calculations  were  performed  with  the  aid  of  Maple  (Appendix  B)  and  a  plot 
of  residuals  in  centimeters  is  shown  in  Figures  (2)  and  (3).  A  linear  change  of  variable  was 
used  to  map  the  trigonometric  polynomial  defined  on  the  interval  [0, 937r/48]  into  the  interval 
1-1.1)- 
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Figure  (2)  -  Trigonometric  polynomial  residuals  (cm) 


Trigonometric  Polynomial  Residuals 
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Figure  (3)  -  Trigonometric  polynomial  residuals  (cm),  (zoomed) 


The  subinterval  chosen  for  this  analysis  was  located  at  the  center  of  the  set  of  nodes  in 
order  to  show  the  location  where  the  trigonometric  interpolant  is  most  accurate.  Examina- 
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tion  of  Figure  (5)  reveals  that  on  a  single  subinterval  the  trigonometric  polynomial  method 
agrees  with  the  Lagrange  interpolation  method  to  within  roughly  ten  meters.  The  Maple 
worksheet  used  to  generate  this  figure  is  included  as  Appendix  C  at  the  end  of  this  paper. 

The  effects  of  error  growth  near  the  ends  of  the  intervals  for  the  FFT  method  could  'be 
handled  by  shifting  the  twenty-four  hour  period  over  which  the  trigonometric  polynomial  is 
derived,  placing  the  data  point  squarely  in  the  center  for  the  most  accurate  work.  A  bound 
on  the  mean  square  error  incurred  by  approximating  the  periodic  function  /  from  which  the 
data  is  sampled  by  the  interpolant  p  is  given  in  Briggs  and  Henson  [1995]  as 

11-^  ~  —  yy-p+i 

where  N  is  the  number  of  data  points,  C  is  &  constant  and  the  periodic  extension  of  /  has 
(p  —  1)  continuous  derivatives,  p  >  1.  Since  continuity  of  the  periodic  extension  of  /  is 
required  for  this  bound  we  cannot  use  it  unless  the  function  is  truly  periodic  on  the  chosen 
interval.  Briggs  and  Henson  [1995]  perform  a  trigonometric  interpolation  on  an  arbitrary 
polynomial  defined  on  [—1, 1]  and  comment  that  it  is  not  unreasonable  to  suspect  the  mean 
square  errors  to  decrease  as  N~^. 

Another  concern  will  be  the  time  required  to  evaluate  the  interpolant,  so  it  is  of  interest 
to  discuss  the  rapid  calculation  of  the  trigonometric  polynomials,  in  particular  the  terms 

cos  TTt,  sin  wt,  cos  27rf,  sin  27rt, . . . ,  cos  nnt,  sin  mrt. 

Using  the  trigonometric  summation  formulas  one  can  write  the  well-known  recursion  [Go- 
ertzel,  1960] 

CTjk  =  cos  kirt 
Tk  =  sin  kirt 

{Z)-in7){n)  '  = 

which  requires  only  two  trigonometric  calculations  in  order  to  recover  all  the  needed  terms. 
■Of  course,  the  growth  of  round-off  errors  should  always  be  checked  when  using  a  long  sequence 
(i.e.  n  large). 

It  is  important  that  we  use  all  the  available  data  in  forming  the  trigonometric  polyno¬ 
mial  pn{t)  because  the  mean  square  error  involved  in  using  the  FFT  method  depends  so 
critically  on  the  number  of  data  points.  Since  96  is  not  a  power  of  two  and  since  the  DFT 
is  computationally  expensive,  we  must  resort  to  another  algorithm. 


Tshebyshev  Polynomial  Interpolation 

In  order  for  us  to  understand  Tshebyshev  polynomial  interpolation,  it  is  necessary  to 
re-examine  the  polynomial  interpolation  of  n  -|- 1  data  points  on  a  given  interval.  Given  n  -f- 1 
points  on  an  interval,  the  nth  degree  polynomial  which  interpolates  those  points  over  the 
given  interval  is  unique.  Of  course  this  does  not  prevent  us  from  writing  the  polynomial  in 
a  number  of  different  ways.  Before  we  describe  what  Tshebyshev  polynomial  interpolation 
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is,  it  may  be  a  good  idea  to  say  what  it  is  not.  Tshebyshev  polynomial  interpolation  is 
not  the  expression  in  the  Tshebyshev  polynomial  basis  of  an  nth  order  polynomial  which 
interpolates  n  + 1  arbitrary  data  points  on  a  given  interval.  We  are  free  to  choose  any  suitable 
polynomial  basis  for  expressing  such  a  polynomial,  but  we  must  remember  that  the  results 
of  any  calculations  performed  using  that  polynomial  will  be  the  same  regardless  of  how  the 
polynomial  is  expressed;  simply  writing  a  given  polynomial  in  a  different  basis  does  not  alter 

it.  Therefore  the  residuals  which  are  produced  when  using  these  different  forms  of  the  same 
polynomial  will  be  the  same.  However,  another  set  of  n  +  1  points  on  the  same  interval 
would  yield  a  different  polynomial.  As  it  turns  out,  polynomial  interpolation  is  sensitive  to 
the  distribution  of  the  points  being  interpolated.  If  we  were  allowed  to  choose  the  points  on 
an  interval  to  be  interpolated,  we  would  find  that  certain  choices  of  n  +  1  points  would  yield 
polynomials  which  did  a  better  job  of  interpolating  than  others.  As  stated,  our  data  points 
are  evenly  spaced  throughout  the  interval  of  interest,  so  that  will  not  have  the  luxury  of 
choosing  a  preferred  set  of  n  +  1  points  on  any  interval  unless  we  are  able  to  generate  more 
data.  The  following  discussion  closely  follows  the  treatment  given  by  Press  et  al  [1992].  The 
Tshebyshev  polynomial  of  degree  n  on  [—1,1]  is  defined  as 

Tn{t)  =  cos(n  arccos  f) 

It  follows  that  the  Tshebyshev  polynomials  satisfy  the  three  term  recurrence  relation 


In  addition, 


and 


/ 


/: 


Tnit)Tkit) 

-1  a/1 

Tnit)Tk{t) 


Vn-lit) 

n  =  1, 2, 

=  0, 

k 

fTT 

A;  —  n  =  0 

U/2 

A;  =  n  ^  0 

so  that  the  Tshebyshev  polynomials  are  orthogonal  on  the  interval  [—1,1]  with  respect  to 
the  weight  function  w(t)  =  l/y/l  —  The  first  few  Tshebyshev  polynomials  are  given  by 

To{t)  =  1 


n{t)  =  t 
T^it)  =  2t^-l 
Tz{t)  =  At^  -  3f 
T^{t)  =  -8t^  +  I 

Tsit)  =  16t®  -20t^  +  5f 


It  can  be  shown  that  the  zeros  of  Tn{t)  on  [—1, 1]  are  given  by 

\{j-l/2y 


tj  =  COS 


n 


j  =  l,2,...,n 


14 


and  that  the  Tshebyshev  polynomials  satisfy  a  discrete  orthogonality  condition  [Press  et  al, 
1992] 

^  0  ^  k 

=  <  n/2  i  =  A;  0 

V  n  i  =  A:  =  0 

The  powers  of  t  can  also  be  expressed  in  the  Tshebyshev  polynomial  basis. 

Any  function  f{t)  can  be  approximated  by  a  finite  linear  combination  of  Tshebyshev 
polynomials,  i.e.  ^ 

f{t)  =  -^0,0  + 

^  »=i 


where 


2  /pO,, 

IT  J-1  v/1  — 


These  coefficients  can  be  computed  numerically  using  the  M  equally  spaced  data  points 

tj  =  -l+jAt,  j  =  ,M  -1 


M-1 

e.g.  via  the  composite  midpoint  rule  (we  have  to  avoid  evaluating  the  integrand  at  the 
endpoints,  to  and  tu-i) 


^  j=o 


^  ^i+i/2 


where  f{tj+i/2)  can  be  approximated  by  a  quadratic  or  cubic  polynomial.  For  example  to 
use  cubic  splines  to  approximate  the  integrand,  we  get 

O  r  A4  M-2  l3  M-2  l2  M-2  M-2 

^  4  o  I  j-o 

where  the  integrand  g{t)  is  approximated  on  the  interval  [tj,tj+i]  by  the  cubic  polynomial 

g{t)  =  aj{t  -  tjf  +  bj{t  -  tjf  +  Cj{t  -  tj)  +  dj 

We  can  also  use  the  midpoint  rule  on  the  two  leftmost  panels  and  the  two  rightmost 
panels  only.  On  the  rest  of  the  panels  we  can  use  the  fourth  order  Simpson’s  4  rule.  This 
will  yield  the  following  approximation  for  the  coefficients: 

AA  r  284  8  ,2  ,  ,  1 

a,-  =  —  \4g2  +  -gs  +  -^g4  +  -^95 - H  -^gM-s  +  ■^gM-2  +  4flfAf-i  J  , 


where 


f{tj)Ti(tj) 


h  - 1] 
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take  the  trapezoidal 

It::; Xlu  :  odd  oun.be.  of  panels  e 

To  evaluate  f{r)  we  then  use  the  following  recurs 

djV+l  ~  ci|v+2 

o  =  iV,N  - 

d^=2rd,+i--dj+2  +  ^v  ^ 

ji^r)  =  do  =  ^di  ^2+2  elation  on  the  Interval  is 


Tlta.. 

“tnhtttrefhods  should  be  — <>• 

run  times  ot  tnese 
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Appendix  A 


This  section  provides  a  FORTRAN  subroutine  which  can  be  used  to  perform  the  inter¬ 
polation  discussed  in  this  paper. 


Program  lagreinge 
c 

c  implements  lagrange  interpolation  of  a  function 

c  given  as  equally  spaced  data  points 

c 

c  let  t.j=-l+j  Delta  t,  Delta  t=2/(M-l),  M=96 

c 

implicit  real*8  (a-h,o-z) 
real+8  t(96) ,f (96) ,y(96) 
real *4  tcLray(2)  ,paray(2) 
atony=dtime (taray) 

c  m  is  number  of  data  points 

m=96 

c  mp  is  the  number  of  interpolating  points 
mp=12 

print  m=  ',m,'  mp=',mp 
c 

c  get  m(=96)  equally  spaced  points  on  [-1,1] 

c 

deltat=2 . / df loat (m-1) 


print  *,  'delta  t  ' ,deltat 

t(l)=-l. 

do  10  j=2,m 

10  t(j)=t(j-l)+deltat 

c 

c  compute  f  values 

c  synthetic  data 

c 

call  f inter(t ,m,f ) 
c 

c  print  the  points 
c 

print  12,(j,t(j),f(j),j=l,m) 
12  format(2x,i3,2el4.7) 

c 
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c  compute  y  values 

c 

call  yi(t,f,mp,y) 
atony=etime(taray) 

write(6,*)  '  User  time  =  taray(l)  seconds’ 
write(6,*)  ’  System  time  =  taray(2) , ’  seconds’ 
write(6,*)  ’  Total  time  =  ’,  taray(l)  +  taray(2),’  seconds’ 
bt  0  ny =dt ime ( par ay ) 

ttt=- 1 . dO- . 5d0*deltat 

print  t  ’.’approximate  exact  error’ 

do  5000  ij=l,mp 
ttt=ttt+deltat 
c 

c  compute  phi_n  (tau) 

c 

call  phi(t,ttt .mp.phin) 
sum=0 . 

do  1000  j=l,mp 

1000  sum*sum+y(j)/(ttt-t(j)) 

pn=phin*sum 
call  finter(ttt,l,ddd) 
print  *,ttt ,pn,ddd,pn-ddd 
5000  continue 

btony=etime(peiray) 

write(6,*)  ’  User  time  =  paray(l)  ,’  seconds’ 
write(6,*)  ’  System  time  =  ’,  paray(2),’  seconds’ 
write(6,*)  '  Total  time  =  ’,  paray(l)  +  paray(2),’  seconds’ 

stop 

end 

subroutine  f inter(t ,m,f) 
implicit  real *8  (a-h,o-z) 
real*8  t(l) ,f (1) 
do  10  j=l,m 

f  (j)=dsin(t(j))+dcos(t(j)) 

10  continue 

return 
end 

subroutine  yi(t,f,m,y) 
implicit  real*8  (a-h,o-z) 
real*8  t(l)  ,f (1) ,y(l) 
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do  20  i=l,m 

y(i)=f(i) 
do  10  j  =  l,in 
if (j . eq. i)  go  to  10 
tt=t(i)-t(j) 
y(i)=y(i)/tt 
10  continue 

20  continue 

return 
end 

subroutine  phi(t ,ttt,m,phin) 
implicit  real*8  (a-h,o-z) 
real*8  t(l) 
phin=l .dO 
do  10  j=l,m 
phin=phin* (ttt-t ( j ) ) 

10  continue 

return 
end 
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Appendix  B 


This  section  gives  a  Maple  program  for  Lagrange/FFT  Comparison  using  GPS  Ephemeris. 
\#  Christopher  P.  Sagovac,  29  May  1995 

\#  The  following  x-coordinate  data  from  satellite  1  in  the  file  jpl07365.sp3 
is  used  to  compare  the  Lagrange  Polynomial  and  FFT  methods  of 
interpolation.  96  points  separated  at  intervals  of  900  seconds  cover  23hrs 
45mins.  11th  order  Lagrange  interpolation  is  performed  on  12  evenly 
distributed  points  on  a  subinterval  within  the  ephemeris;  remaining  points 
are  used  to  generate  residuals.  The  user  specifies  the  subinterval  over 
which  the  interpolation  takes  place  by  setting  the  starting  node  number  in 
the  range  a=1..61.  For  the  Lagrange  interpolation,  the  ephemeris  data  is 
mapped  into  the  interval  [-1,1].  Selected  points  aire  separated  by  2700 
seconds.  The  FFT  is  performed  over  the  entire  ephemeris  using  32  equally 
spaced  points  with  subinterval  spacing  of  2700  seconds;  remaining  points 
are  used  to  generate  residuals.  For  the  FFT  interpolation,  the  ephemeris 
data  is  mapped  into  the  interval  [0,95*$\pi$/48] .  This  allows  a  comparison 
of  the  Lagrange  interpolating  polynomial  and  the  FFT  method  over  an  ephemeris 
interval . 

Residual  comparison  takes  place  at  corresponding  node  numbers  on  each  of  the 
different  sub int erval s . 

First  we  generate  a  list  of  the  ninety-six  abscissas  for  both  intervals: 

>  listl :=[seq(m*2*$\pi$/96,  m=0..95)]: 

>  list2:=[seq(-l+m*2/95,  m=0..95)]: 

\#  A  list  of  the  corresponding  ninety-six  ordinate  values  tciken  from  file 
jpl07365 . sp3 

\#  (these  are  the  same  for  both  sets  of  abscissas) : 

>  satlist2:= 

[15500 .414560 , 15268 . 853410 , 14843 . 274690 , 14200 . 160540 , 13322 .919650 , 

> 12202 . 689520 , 10838 . 841980 , 9239 . 173760 , 7419 . 776750 , 5404 . 594490 , 3224 . 683250 , 
>917.206990,-1475.795270,-3908.822870,-6334.455100,-8704.963870,-10973.938070, 
>-13097 . 860380,-15037 . 578740,-16759 .618700,-18237 . 288490 , -19451 . 535890 , 
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>-20391 . 524780 , -21054 . 909440 , -21447 . 795240 , -21584 .386110, -21486 . 330880 , 

>-21181 . 791750 , -20704.269190 , -20091 . 227140 , -19382 . 570580 , -18619 . 033950 , 

>-17840 . 542760 , -17084 .613010 , -16384 . 851470 , -15769 . 617150 , -15260 . 897820 , 

>-14873 . 447530 , -14614 . 220590,-14482 . 125710 , -14468 .110980 , -14555 . 576720 , 

>-14721 . 100000 , - 14935 . 441300 , -15164 . 792660 , -15372 . 216220 , -15519 . 214910 , 

>-15567 . 371330 , -15479 . 988750 , -15223 . 668540 , -14769 . 761120 , -14095 . 633690 , 

>-13185 . 705450 , -12032 .211050,-1063 . 663960 , -9005 . 003600 , -7157.422240 , 

>-5117. 879990 ,-2918. 327780 , -596 . 668940 , 1804 . 500890 , 4239 . 327820 , 6660 . 254360 , 

>9019 . 637090 ,11271. 368100 , 13372.440740 , 15284 . 402460 , 16974 . 641660 , 

>18417 .461210, 19594 . 898880 , 20497 . 263810 ,21123. 368570 , 21480 . 447000 , 

>21583 . 759770 , 21455 . 901300 , 21125 . 832970 , 20627 . 678050 , 19999 . 323680 , 

>19280 . 882770 , 18513 . 075050 , 17735 . 590120 , 16985 . 496810 , 16295 . 762020 , 

>15693 . 938430 , 15201 . 073980 , 14830 . 887900 , 14589 . 247140 , 14473 . 965140 , 

> 14474 . 931940 , 14574 . 570650 , 14748 . 602380 , 14967 . 088360 , 15195 . 707160 , 

>15397 . 214850 , 15533 . 028660 , 15564 . 870070] : 

\#  Normalize  the  data  by  dividing  by  $10“5$: 

>  list3:=[seq(evalf (satlist2[m]/100000,ll) ,m=l. .96)] : 

\#  We  next  form  the  11th  order  Lagrange  interpolating  polynomial  using  twelve 
evenly  spaced  points  on  the  interval.  The  other  unused  points  will  be  used  to 
examine  residuals.  Setting  the  parameter  "a"  (below)  chooses  the  point  at  which 
the  interval  begins;  valid  choices  axe  a=1..61.  Subintervals  chosen  to  lie 
neax  the  ends  of  the  ephemeris  interval  represent  a  worst  case  for  the  FFT 
method,  which  attains  its  greatest  accuracy  in  the  center  of  the  interval. 

For  instance,  the  leftmost  interval  of  length  36  points  corresponds  to  a=l. 

The  rightmost  interval  of  length  36  corresponds  to  a=61.  In  this  way  one  can 
peruse  the  entire  ephemeris  interval  to  see  how  the  32  point  FFT  method 
compaxes  to  the  local  step  along  Lagrange  method.  Choose  a  value  in  the 
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center  of  the  interval,  a=32: 


>  a: =32: 

1 

\#  Here  are  the  values  of  the  abscissas  for  the  chosen  subinterval  on  the 
interval  [-1,1): 

>  plistl :=[seq(evalf (list2[a+3*m] , 16) ,  m=0..11)]: 

\#  Here  are  the  values  of  the  ordinates  (nodes)  for  the  chosen  subinterval 
(from  list3)  : 

>  plist2:=Cseq(list3[a+3*m] ,m=0. . 11)] : 

\#  Now  calculate  the  11th  degree  Lagrange  interpolating  polynomial  over  the 
chosen  subinterval: 

>  p:=interp(plistl,  plist2,  z) : 

>  readlib (FFT) : 

\#  Form  a  list  of  32  evenly  distributed  data  points  from  list3: 

>  fftlist2:  =  [seq(list3[l+3*m] ,m=0. .31)]  : 

>  X  :=  array (fftlist2) : 

\#  There  is  no  imaginary  component  for  this  real  data  set;  y  is  a  zero  array 


>  y  :  = 

array ([0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0,0, 

\#  Use  hardware  floating  point  number  system  for  fast  calculation: 

>  evalhf (FFT(5,var(x) ,var(y))) : 

>  r :=map(x->x/32, convert (x,  list)): 

>  s:=map(y->y/32, convert (y, list)) : 

>  evenf  :=r[l]+r[17]*cos(16*z) : 

>  for  n  from  1  to  15  do  evenf :=  evenf+2*r[n+l]*cos(n*z)  od: 


0]): 
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>  oddf:=0: 


>  for  n  from  1  to  15  do  oddf:=oddf+  (-2) *s [n+l] *sin(n*z)  od: 

\#  Here  then  is  the  interpolating  trigonometric  polynomial,  good  for  the 
entire  24  hour  period: 

>  g : =oddf +evenf : 

\#  Plot  the  difference  between  the  FFT  method  and  the  Lagrange  interpolating 
polynomial  on  [-1,1]: 

>  FFTResidsl :=[seq( [-.3+m/150,$10"{10}$*abs(evalf (subs(z=evalf (- .3+m/l50 ,20) , 

p),16)- 

>  evalf (subs(z=95/96*$\pi$*( .7+m/150) ,  g) ,20))] ,m=0. .90)] : 

\#  interf ace(plotdevice=postscript ,  plotoutput=interp4) ; 

>  plot (FFTResidsl ,  style=point,  title= ‘Trigonometric  Polynomial  Residuals'); 

> 

FFTResids2 : = [seq( [evalf (m/1000) , $10“{103-$*abs (evalf (subs (z=evalf (evalf (m/ 1000) 

,20),  p),16)- 

>  evalf (subs(z  =95/96*$\pi$*(l+m/1000)  ,  g) ,20)3] ,m=0 . .90)] : 

\#  interface (plotdevice=postscript ,  plotoutput=interp5) ; 

>  plot(FFTResids2,  style=point,  title='Trigonometric  Polynomial  Residuals'); 
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Appendix  C 


This  section  gives  a  Fortran  program  for  Tshebyshev  approximation  using  Simpson ’n  | 
rule  for  the  numerical  approximations  of  the  coefficients. 

subroutine  cheb(rr ,a,ti ,n,alf ,t ,f ,m) 
implicit  real*8  (a-h,o-z) 
realms  a(l) ,ti(l) ,alf (1) ,t(l) .f (1) 
c 

c  rr  is  delta  t  /  pi 

c  (t,f)  are  the  equally  spaced  data  points 

c  n  is  the  highest  degree  of  Tshebyshev  polynomial  in  the  expansion 
c  m  is  the  number  of  data  points 
c 

c  a(i)  are  the  coefficients  of  the  expansion  in  Tshebyshev  polynomials 
c 

do  20  i=l,n 
20  a(i)=0. 

c 

c  weights  (alf)  for  midpoint/trapezoidal/simpson  combination 

c 

alf(2)=4.d0 
alf (3)=2.d0/3.d0 
do  22  i=4,m-4,2 
alf(i)=8.d0/3.d0 
alf (i+l)=4.d0/3.d0 
22  continue 

alf (m-3) =5 . dO/3 . dO 
alf (m-2)=l.d0 
alf (m-l)=4.d0 
c 

c  integration 

c 

do  30  j=2,m-2 
tt=t(j) 

den=dsqrt(l .-tt*tt) 

c  print  ’ j ,tt ,den=' , j ,tt ,den 

do  40  i*l,n 

c 

c  compute  all  Tshebyshev  polynomial  up  to  degree  n  at  point  tt 
c 

call  cheby(ti,n,tt) 
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t  enii=t  i(i)*f(j)*alf(j)/den*rr 
a(i)=a(i)+term 

c  print  i  ti(j) ,a(i) ' , j ,i,ti(i) ,a(i) 

40  continue 

30  continue 

c  print  ’  a(i)=  ’  ,  a 

return 

end 
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